High-radix multiplier-divider

ABSTRACT

The high-radix multiplier-divider provides a system and method utilizing an SRT digit recurrence algorithm that provides for simultaneous multiplication and division using a single recurrence relation. When A, B, D and Q are fractions (e.g., Q=0·q −1  q −2  . . . q −n ), then the algorithm provides for computing 
             S   =     AB   D           
to yield a w-bit quotient Q and w-bit remainder R by: (1) determining the next quotient digit q −j  using a quotient digit selection function; (2) generating the product q −j D; and (3) performing the triple addition of rR j-1 , (−q −j D) and
 
               b     -     (     j   -   1     )         ⁡     (     A   r     )           
where R 0 =b −1 A r   −1 . The recurrence relation may be implemented with carry-save adders for computation using bitwise logical operators (AND, OR, XOR).

CROSS-REFERENCE TO RELATED APPLICATION

This application is a continuation-in-part of U.S. patent application Ser. No. 11/819,749, filed Jun. 28, 2007.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to high performance digital arithmetic circuitry, and more particularly to a high-radix multiplier-divider implemented in hardware for efficiently computing combined multiplication and division operations, e.g., computing ((A·B)÷D), as well as modulo multiplication, e.g., (A·B) mod D.

2. Description of the Related Art

Many of the recent computer intensive applications, e.g., multimedia and public-key cryptosystems, have been witnessing a continuous increase in demand for more computational power. For example, public-key cryptosystems make heavy use of the modulo multiplication operation, which comprises a multiplication operation together with a division/reduction operation. The key size of RSA public-key cryptosystems has been continuously getting larger, from 512 bits to 1024 bits, and most recently to 2048 bits, causing increased demand for more computational power. There exists a quite extensive literature that describes the theory and design of high speed multiplication and division algorithms. Division algorithms are divided into five classes: (1) digit recurrence; (2) functional iteration; (3) very high radix; (4) table lookup; and (5) variable latency.

Digit recurrence is the oldest class of high speed division algorithms, and, as a result, a significant quantity of literature has been written proposing digit recurrence algorithms, implementations, and techniques. The most common implementation of digit recurrence division in modern processors has been the Sweeney, Robertson, Tocher (SRT) method.

Digit recurrence algorithms retire a fixed number of quotient bits in every iteration. Implementations of digit recurrence algorithms are typically low complexity, utilize small area, and have relatively large latencies. The fundamental choices in the design of a digit recurrence divider are the radix, the allowed quotient digits, and the representation of the partial remainder (residue). The radix determines how many bits of quotient are retired in an iteration, which fixes the division latency. Larger radices can reduce the latency, but increase the time for each iteration. Judicious choice of the allowed quotient digits can reduce the time for each iteration, but with a corresponding increase in complexity and hardware. Similarly, different representations of the partial remainder (residue) can reduce iteration time, but with corresponding increases in complexity.

Digit recurrence division algorithms use iterative methods to calculate quotients one digit per iteration. SRT is the most common digit recurrence division algorithm. The input operands are represented in a normalized floating point format with w-bit significands in sign and magnitude representation. Assuming floating point representation, the algorithm is applied only to the magnitudes of the significands of the input operands. Techniques for computing the resulting exponent and sign are straightforward. The most common format found in modern computers is the IEEE 754 standard for binary floating point arithmetic. This standard defines single and double precision formats.

In a division operation (N/D), N is a 2 w-bit dividend, while D is a w-bit divisor. The division result is a w-bit quotient Q (Q=0·q⁻¹ q⁻² . . . q_(−n)) and a w-bit remainder R such that N=QD+R and |R|<|D|. The w-bit quotient is defined to consist of n radix-r digits with r=2^(m), (w=n×m). A division algorithm that retires m quotient bits per iteration is said to be a radix-r algorithm. Such an algorithm requires n iterations to compute the result. For no overflow, i.e., so that Q is w-bits, the condition |N|<|D| must be satisfied when dividing fractions. The following recurrence is used in every iteration of the SRT algorithm; R _(j) =rR _(j-1) −q _(−j) D j=1,2,3, . . . , n; where

-   -   R₀=N     -   D=divisor,     -   N=dividend, and     -   R_(j) is the partial residue at the j^(th) iteration.         One quotient digit (m-bits) is retired each iteration using a         quotient digit selection function SEL where:         q_(−j)=SEL(rR_(j-1), D).

After n iterations, the final value of the quotient Q and the remainder R are computed from R_(n) as follows:

${{{if}\mspace{14mu} R_{n}} > 0},{{{then}\mspace{14mu} Q} = {{\sum\limits_{j = 1}^{n}\;{q_{- j}r^{- j}\mspace{14mu}{and}\mspace{14mu} R}} = {R_{n}r^{- n}}}}$ ${{else}\mspace{14mu}\left( {Q = {\sum\limits_{j = 1}^{n}{q_{- j}r^{- j}}}} \right)} - {r^{- n}\mspace{14mu}{where}\mspace{14mu} r^{- n}\mspace{14mu}{is}\mspace{14mu} 1\mspace{14mu}{in}\mspace{14mu}{the}\mspace{14mu}{least}}$ significant  position, or $Q = {\left( {\sum\limits_{j = 1}^{n}{q_{- j}r^{- j}}} \right) - {{ulp}\mspace{14mu}{where}\mspace{14mu}{ulp}\mspace{14mu}{designates}\mspace{14mu} a\mspace{14mu}{unit}\mspace{14mu}{in}\mspace{14mu}{the}}}$ least  significant  position  and R = (R_(n) + D)r^(−n).

The critical path of the basic SRT digit recurrence algorithm comprises the following steps: (1) determination of the next quotient digit q_(−j) using the quotient digit selection function, a look-up table typically implemented as a PLA, or read only memory (ROM); (2) generation of the product q_(−j)D; and (3) subtraction of Q_(−j)D from the shifted partial residue rR_(j-1). Each of these steps contributes to the algorithm cost and performance.

A common method to decrease the overall latency of the algorithm (in machine cycles) is to increase the radix r of the algorithm. Assuming the same quotient precision, the number of iterations of the algorithm required to compute the quotient is reduced by a factor f when the radix is increased from r=2^(m) to r=2^(mf). For example, a radix-4 algorithm retires two bits of quotient in every iteration. Increasing to a radix-16 algorithm allows for retiring four bits in every iteration, halving the latency.

This reduction does not come for free. As the radix increases, the quotient-digit selection becomes more complicated and, accordingly, slower to compute. Since the quotient selection logic is on the critical path of the basic algorithm, using higher radices causes the total cycle time of the division iteration to increase. The number of cycles, however, is reduced for higher radices. As a result, the total time required to compute a w-bit quotient may not be reduced as expected. Furthermore, the generation of divisor multiples may become impractical or infeasible for higher radices. Thus, these two factors can offset some of the performance gained by using higher radices.

Typically, for a system with radix r, a redundant signed digit set (D_(a)) is used to increase the performance of the algorithm. To be redundant, the size of the digit set should be greater than r, including both negative and positive digits. Thus, q_(−j)εD_(α)={−α, −α+1, . . . , −1, 0, 1, . . . , β−1, β}, where the number of allowed digits (α+β+1) is greater than r. It is fairly common to choose a symmetric digit set where α=β, in which case the size of the digit set (2α+1)>r, which implies that a must satisfy the condition α≧┌r/2┐. The degree of redundancy is measured by the value of the redundancy factor h, where h=α/r⁻¹. Redundancy is maximal when α=r−1, in which case h=1, while it is minimal when α=r/2 (i.e., ½<h≦1).

For the computed R_(j) value to be bounded, the value of the quotient digit must be selected such that |R_(j)|<hD. Using larger values of h (i.e., large α) reduces the complexity and latency of the quotient digit selection function. This, however, results in a more complex generation of the divisor multiples. Divisor multiples that are powers of two can be formed by simple shifting, while those that are not powers of two (e.g., three) require additional add/subtract steps. The complexity of the quotient digit selection function and that of the generating divisor multiples must be balanced.

To define the quotient digit selection function, a containment condition is used to determine the selection intervals. A selection interval is the region defined by the values of the shifted partial residue (rR_(j-1)) values and the divisor (D) in which a particular quotient digit may be selected. The selection interval is defined by the upper (U_(k)) and lower (L_(k)) bounds for the shifted partial residue (rR_(j-1)) values in which a value of quotient digit q_(j)=k may be selected to keep the partial residue R_(j) bounded. These are given by: U _(k)=(h+k)D and L _(k)=(−h+k)D.

The P-D diagram is a useful tool in defining the quotient-digit selection function. It plots the shifted partial residue (P=rR_(j-1)) versus the divisor D. The U_(k) and L_(k) straight lines are drawn on this plot to define selection interval bounds for various values of k. FIG. 6 shows a P-D diagram for the case where r=4 and α=2(h=⅔). The shaded regions are the overlap regions where more than one quotient digit may be selected. Table I shows representative values of the upper and lower bounds U_(k) and L_(k) for the permissible values of the quotient digit k.

TABLE I Upper and Lower Bounds vs. Quotient Digit k U_(k) = (h + k)D L_(k) = (−h + k)D −2 −4/3D   −8/3D −1 −1/3D   −5/3D 0 2/3D −2/3D 1 5/3D   1/3D 2 8/3D   4/3D

There is a need for a digital multiplier-divider unit that can efficiently compute

$S = \left( \frac{A \cdot B}{D} \right)$ where the multiplicand A, the multiplier B, and the divisor D are w-bit unsigned numbers. Computing S yields a w-bit quotient and a w-bit remainder R such that: A·B=Q·D+R and |R|<|D|.

Conventionally, S would be computed using two independent operations: a multiplication operation, and a division operation. Whereas digit recurrence relations for these two operations have been proposed and are in common use by digital processors, no single recurrence relation has been proposed to simultaneously perform the multiplication and division operations as needed to efficiently compute

$S = {\left( \frac{AB}{D} \right).}$

Thus, a high-radix multiplier-divider solving the aforementioned problems is desired.

SUMMARY OF THE INVENTION

The high-radix multiplier-divider provides a system and method utilizing an SRT digit recurrence algorithm that provides for simultaneous multiplication and division using a single recurrence relation. When A, B, D and Q are fractions (e.g., Q=0·q⁻¹ q⁻² . . . q_(−n)), then the algorithm provides for computing

$S = \frac{AB}{D}$ to yield a w-bit quotient Q and w-bit remainder R by: (1) determining the next quotient digit q_(−j) using a quotient digit selection function; (2) generating the product q_(−j)D; and (3) performing the triple addition of rR_(j-1), (−q_(−j)D) and

$b_{- {({j - 1})}}\left( \frac{A}{r} \right)$ where R₀=b⁻¹Ar⁻¹. The recurrence relation may be implemented with carry-save adders for computation using bitwise logical operators (AND, OR, XOR).

These and other features of the present invention will become readily apparent upon further review of the following specification and drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of an exemplary circuit implementing a high-radix multiplier-divider according to the present invention.

FIG. 2 is chart showing an exemplary multiplier-divider operation according to the present invention.

FIG. 3 is a comparison constants graph defining quotient digit values in a high-radix multiplier-divider according to the present invention.

FIG. 4 is a P-D diagram showing overlap regions for P>0 in a high-radix multiplier-divider according to the present invention.

FIG. 5 is a P-D diagram showing overlap regions for P<0 in a high-radix multiplier-divider according to the present invention.

FIG. 6 is a P-D diagram for P>0 with specified r and alpha values according to the prior art.

FIG. 7 is a table showing a set of solution steps in a high-radix multiplier-divider according to the present invention.

FIG. 8 is a continuation of the table of FIG. 7 showing the set of solution steps in a high-radix multiplier-divider according to the present invention.

FIG. 9 is a block diagram illustrating a computer system for integrating and performing the high-radix multiplier-divider.

Similar reference characters denote corresponding features consistently throughout the attached drawings.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention provides a digit recurrence (SRT) multiplier-divider apparatus and method that utilizes a single recurrence relation. The method includes an algorithm that may be implemented in software, but is preferably implemented in hardware for greater speed.

The apparatus includes a circuit configured to carry out the algorithm. The circuit may be incorporated into the architecture of a computer processor, into a security coprocessor integrated on a motherboard with a main microprocessor, into a digital signal processor, into an application specific integrated circuit (ASIC), or other circuitry associated with a computer, electronic calculator, or the like. The method may be modified so that the circuit may include carry-propagate adders, or the circuit may include carry-save adders, or it may include compressors, e.g., (4-2) compressors.

The method can perform simultaneous multiplication and division. Roth the apparatus and the method may be utilized in a variety of applications, including but not limited to, networked computers, public-key cryptosystems, multimedia applications, digital communication devices, and the like, where the method and circuitry provide for high speed performance of modular arithmetic operations involved in the encryption and decryption of messages, where the method and the circuitry provide increased speed for greater circuit efficiency, increased productivity, and lower network, processor, system overload and costs.

It is desired to provide a digital multiplier-divider unit that can efficiently compute

$S = \left( \frac{A \cdot B}{D} \right)$ where the multiplicand A, the multiplier B, and the divisor D are w-bit unsigned numbers. Computing S yields a w-bit quotient Q and a w-bit remainder R such that: A·B=Q·D+R;  (1) |R|<|D|.  (2)

In one embodiment, in order to speed-up the computation

${S = \left( \frac{A \cdot B}{D} \right)},$ the recurrence relation of the present invention uses a high radix r=2^(m) where m≧1. The operands A, B, and D are n-digit integers, i.e., A=(a_(n−1), . . . a₁, a₀), B=(b_(n−1), . . . b₁, b₀), and D=(d_(n−1), . . . d₁, d₀), where

$n = \left\lceil \frac{w}{m} \right\rceil$ and a_(i), b_(i), d_(i)ε{0, 1, . . . , r−1}.

Thus, the present invention provides for an enhanced multiply-divide recurrence relation given by: R _(j) =rR _(j-1) −q _(n−j) Dr ^(n) +b _(n−j-1) Ar ^(n−1) j=1, 2, . . . , n  (3)

where,

q_(i) is the i^(th) quotient digit

b_(i) is the i^(th) digit of B

b⁻¹=0,

R_(j) is the j^(th) running partial remainder and

R₀=b_(n−1)Ar^(n−1).

The final quotient Q and remainder R results are given by:

$\begin{matrix} {{Q = {q_{n - 1}q_{n - 2}\mspace{14mu}\ldots\mspace{14mu} q_{2}q_{1}q_{0}}},{i.e.},{Q = {\sum\limits_{j = 0}^{n - 1}\;{q_{j}{r^{j}.}}}}} & (4) \end{matrix}$ If R_(n)<0 then the following correction step is performed:

Q=Q−ulp, where ulp designates a unit in the least significant position,

and R_(n)=R_(n)+D, with

$\begin{matrix} {R = \frac{R_{n}}{r^{n}}} & (5) \end{matrix}$

The following shows that executing the n iterations of the proposed recurrence relation yields the desired Q & R values as defined by equations (1) and (2):

  R_(j) = rR_(j − 1) − q_(n − j)Dr^(n) + b_(n − j − 1)Ar^(n − 1)   R₀ = b_(n − 1)Ar^(n − 1) $\mspace{20mu}\begin{matrix} {R_{1} = {{r^{n}b_{n - 1}A} - {q_{n - 1}{Dr}^{n}} + {r^{n - 1}b_{n - 2}A}}} \\ {= {{A\left( {{r^{n}b_{n - 1}} + {r^{n - 1}b_{n - 2}}} \right)} - {{Dr}^{n}q_{n - 1}}}} \end{matrix}$ $\mspace{20mu}\begin{matrix} {R_{2} = {{A\left( {{r^{n + 1}b_{n - 1}} + {r^{n}b_{n - 2}}} \right)} - {{Dr}^{n + 1}q_{n - 1}} - {q_{n - 2}{Dr}^{n}} + {b_{n - 3}{Ar}^{n - 1}}}} \\ {= {{A\left( {{r^{n + 1}b_{n - 1}} + {r^{n}b_{n - 2}} + {r^{n - 1}b_{n - 3}}} \right)} - {D\left( {{r^{n + 1}q_{n - 1}} + {r^{n}q_{n - 2}}} \right)}}} \end{matrix}$ R₃ = A(r^(n + 2)b_(n − 1) + r^(n + 1)b_(n − 2) + r^(n)b_(n − 3) + r^(n − 1)b_(n − 4)) − D(r^(n + 2)q_(n − 1) + r^(n + 1)q_(n − 2) + r^(n)q_(n − 3))   ⋮ $\mspace{20mu}{{R_{n} = {{{A{\sum\limits_{i = 1}^{n + 1}\;{r^{{2\; n} - i}b_{n - i}}}} - {D{\sum\limits_{j = 1}^{n}\;{r^{{2\; n} - j}q_{n - j}\mspace{14mu}{with}\mspace{14mu} b_{- 1}}}}} = 0}},\mspace{20mu}{R_{n} = {{r^{n}\left( {{AB} - {DQ}} \right)}.}}}$

Accordingly,

$R = {\frac{R_{n}}{r^{n}} = {{AB} - {DQ}}}$ and accordingly, AB=DQ+R.

If the digits of Q are chosen such that the magnitude of the partial residue R_(j) is maintained less than the magnitude of D, then Q is effectively the required quotient of the division operation: AB/D. Since AB=DQ+R and |R|<|D|, then R is indeed the division operation final remainder.

The previous recurrence relation (equation 3) can be rewritten assuming A, B, D, and Q to be fractions of the form B=0·b⁻¹ b⁻² . . . b_(−n), and Q=0·q⁻¹ q⁻² . . . q_(−n). This does not change the computation procedure in any way, and integer operations can be readily mapped to the fractional form. The fractional formulas are more convenient in mathematical representation, however, since they are readily adaptable to floating point representations. The fractional form is obtained from the integer form as follows: A=A _(integer) *r ^(−n) B=B _(integer) *r ^(−n) D=D _(integer) *r ^(−n) R=R _(integer) *r ^(−n).  (6)

Following is the modified fractional multiply-divide recurrence relation:

$\begin{matrix} {{{R_{j} = {{{rR}_{j - 1} - {q_{- j}D} + {b_{{- j} - 1}{Ar}^{- 1}\mspace{31mu}{for}\mspace{14mu} j}} = 1}},2,\ldots\mspace{14mu},n}{{where},{R_{0} = {b_{- 1}{Ar}^{- 1}}}}} & (7) \\ {b_{- i} = {{0\mspace{14mu}{for}\mspace{14mu} i} > n}} & (8) \\ {{Q = {0.q_{- 1}q_{- 2}\mspace{14mu}\ldots\mspace{14mu} q_{- n}}},{i.e.},{Q = {\sum\limits_{j = 1}^{n}{q_{- j}r^{- j}}}}} & (9) \\ {R = {\frac{R_{n}}{r^{n}}.}} & (10) \end{matrix}$

Alternatively, the same recurrence relation may be used with R₀=0. In this case, an extra iteration step of the recurrence relation is needed. Thus; R _(j) =rR _(j-1) −q _(−j) D+b _(−j−1) Ar ⁻¹ for j=1, 2, . . . , n+1  (11) R ₀=0 q ₀=0 b _(−i)=0 for i>n R _(n+1) =r ^(n) R.  (12) The final quotient Q and remainder R are given by:

$\begin{matrix} {{Q = {0.q_{- 1}q_{- 2}\mspace{14mu}\ldots\mspace{14mu} q_{- n}}},{i.e.},{Q = {\sum\limits_{j = 1}^{n}{q_{- j}r^{- j}}}}} & (13) \\ {R = {\frac{R_{n + 1}}{r^{n}}.}} & (14) \end{matrix}$

The previous formulae can be implemented in hardware using shift and add operations. Although the operand size is w-bits (w=mn), the minimum possible size in radix r (r=2^(m)) implementation is w+2 m+1; w-bits for the actual operand, m-bits for the left shift required by the first term (rR_(j-1)), another m-bits for the right shift required by the third term (b_(−j)Ar⁻¹), and finally one extra bit for the sign. FIG. 2 shows an 8-bit binary example 200 (r=2, m=1, and w=n=8) of the multiplier-divider. This example follows exactly the recurrence relation, with the “(shifted)” lines indicating the shifted partial remainders (rR_(j-1)). Note that the integer operands were mapped to their fractional form by simply multiplying each input operand by 2^(−w). Hereinafter, fractional operands are assumed, with the understanding that integer operations can be directly mapped to the fractional form. The width of the operation is 8+2*1+1=11 bits. The quotient digits (bits in this case) were chosen such that the remainder will always lie in the range (−D, +D) and not just [0, +D). Allowing negative remainders can cause the selected quotient digit to be negative. Additionally, in the example of FIG. 2, since all operands (A, B and D) are positive, the final remainder should be positive as well. Since we are allowing negative remainders (with magnitude less than D), we may require a final correction step if the final remainder turns out to be negative. According to the present invention, the correction step adds the divisor D to the partial remainder R_(n) with a corresponding correction to the quotient value by subtracting 1 in the least significant bit position, which is typically referred to as “unit in the least position”, or ulp.

Just as in the case of division, we must have AB<D to guarantee that no overflow may occur (since AB=DQ+R).

The following analysis assumes the use of equations 7-10 with the understanding that similar analysis holds true if the alternative formulae of equations 11-14 are used instead.

Referring to the recurrence relation of equation (7), each iteration includes the following steps: (1) determining the next quotient digit q_(−j) using some quotient digit selection function (q_(−j)=SEL(rR_(j-1), D), which may typically be implemented as a look-up table or as a PLA; (2) generating the product q_(−j)D; and (3) performing the triple addition of rR_(j-1), (−q_(−j)D) and

${b_{- {({j + 1})}}\left( \frac{A}{r} \right)}.$ The resulting partial residue (R_(j)) must guarantee that |R_(j)|<|D|. Satisfaction of the condition |R_(j)|<|D| depends on the proper choice of the quotient digit (q_(−j)).

When performing a multiply-divide operation, we are adding a multiple of the input operand A in each step. The resulting residue thus obtained (R_(j)) cannot be known as predictably as in the case of high radix division. However, we may still restrict the value range of (R_(j)) by placing some restrictions on the value of A. One possible restriction is to impose the constraint that |A|<|D|. This restricts the range of R_(j) as follows: Let A=ωD where, ω<1. R _(j) =rR _(j-1) −q _(−j) D+b _(−j−1) A ^(r−1)

Since Max(b_(i))=(r−1), we have:

${{R_{jmax}\left( q_{- j} \right)} = {{rR}_{j - 1} - {q_{- j}D} + {\left( \frac{r - 1}{r} \right)\omega\; D}}},{or}$ ${{R_{jmax}\left( q_{- j} \right)} = {{rR}_{j - 1} - {q_{- j}D} + {\omega\;{TD}\mspace{14mu}{where}}}},{T = {\left( \frac{r - 1}{r} \right) < 1}},{{{{and}\mspace{14mu}\omega} < 1}\mspace{20mu} = {R_{j{({division})}} + {\omega\;{TD}}}},$ where R_(j(division)) is the residue of a regular high radix division. This shows that the deviation in the remainder curve of the Robertson diagram from the case of pure division can be as high as

$\left( \frac{r - 1}{r} \right)\omega\;{D.}$

The minimum value of R_(j), however, is independent of A and is given by: R _(j min)(q _(−j))=rR _(j-1) −q _(−j) D=R _(j(division)).

Since

$\omega = \left( \frac{A}{D} \right)$ its upper bound value given by A_(max)/D_(min) must be less than 1.

To guarantee satisfaction of this constraint, a pre-processing step shifting A by z-bits to the right is performed. Thus, if the input operand is A′ processing is actually performed on A=A′/2^(z) rather than the input operand A′ itself. Accordingly, the method of the present invention computes S=AB/D and a post-processing step will finally compute S′=A′B/D=S*A′/A=S2^(z). For floating point number representation, an x-bit normalized signific and will result in a ratio of A_(max)/D_(min) that equals [2−2^(−(x−1))] and accordingly, the upper bound of ω is given by:

${\omega \leq \frac{A_{\max}/D_{\min}}{2^{z}}} = {\frac{2\left( {1 - 2^{- {nm}}} \right)}{2^{z}} < {2^{1 - z}.}}$

Since for typical values of n and m, the quantity (1−2^(−nm))<<1, we define the parameter ω_(max)=2^(1-z) as the upper bound for ω such that ω<ω_(max) where: ω_(max)=2^(1-z)≦1.  (15)

The multiply-divide recurrence relation can be implemented in hardware using shift and add operations. Although the problem size is w bits (w=nm), the minimum possible size in radix-r implementations is [(n+2)m+z+1] bits where r=2^(m). Referring to the high-radix multiplier-divider recurrence relation (equation 7), a total of n-digits are needed to accommodate the input operand size, two more digits are needed to account for the left and right shifts (rR_(j-1) & b_(−j−1)Ar⁻¹), z extra bits are needed since computations are performed on the constrained parameter A (A=A′/2^(z)) rather than the input operand A′, and a sign bit is required, since the partial residue R_(j) may be either positive or negative.

Multiplication-division can alternately be performed using a dedicated w×w multiplier producing the 2 w-bit product A×B followed by a dedicated divider to divide the resulting product by D. In addition to the dedicated w×w multiplier, the divider hardware requires adders of [2w+m+1] bits. As an example, consider the case of w=52 bits and r=16 (i.e., m=4 and n=13), with a value of z=8 the multiplier-divider of the present invention requires adders of only 69-bits. The alternate solution of using a dedicated multiplier followed by dedicated conventional divider hardware requires both a 52×52 multiplier, and a divider which uses adders of 109-bits. Furthermore, with floating point input operands, the merged operation of multiplication-division of the present invention requires only one rounding operation. The alternate solution with a dedicated multiplication operation followed by a division operation requires two rounding operations.

Due to the pre-processing step where the input operand A′ is shifted right by Z bit positions, i.e., A=A′/2^(z), a post-processing step where the result S is shifted left by Z bit positions is needed, i.e., S′=S2^(z). In other words, since the resulting quotient and remainder values (Q & R) satisfy the relation AB=QD+R, i.e., (A′*2^(−z))B=QD+R, the true quotient Q′ and remainder R′ that satisfy A′B=Q′D+R′ are computed in a post-processing step by: Q′=Q*2^(z), and R=R*2^(z).

Thus, it is expected that the first Z bits of the resulting quotient (Q) will be zeros. Accordingly, if n-significant digits of Q′ are required, the number of required iterations of the recurrence relation (equation 7) must be raised to

$n + {\left\lceil \frac{Z}{m} \right\rceil.}$

To define the quotient digit selection function, we need to determine the upper and lower bounds of the shifted partial residue (P=rR_(j-1)) for which a given quotient digit value may be selected such that |R_(j)|<|D|. The assumptions under which these bounds can be derived are:

-   -   1. R_(j) is kept bounded, i.e., |R_(j)|<|D|, by defining the         negative and positive range limiting factors h⁻ and h⁺ such         that:         −h ⁻ D≦R _(j) ≦+h ⁺ D where h ⁻ , h ⁺<1.     -   2. The radix r is a power of 2, i.e., r=2^(m).     -   3. The magnitude of A is smaller than the magnitude of D (A=ωD),         where ω<ω_(max)=2^(1-z)≦1).     -   4. The operand B is represented in nonredundant binary format         with radix r digits, i.e., 0≦b_(i)≦r−1. The analysis can be         easily extended to accommodate signed-digit representations of         B.     -   5. A balanced signed redundant quotient digit set D_(q) is used         where:         D _(q)={−α,−α+1, . . . , −1,0,1, . . . ,α−1,α},         with         r/2≦α≦(r−1).     -   6. For the j^(th) iteration, a valid choice of g_(−j) is one         which satisfies the condition −h⁻D≦R_(j)≦+h⁺D where h⁻, h⁺<1.

For a feasible implementation of the high-radix multiplier-divider recurrence relation (equation 7), when the shifted partial residue rR_(j-1) equals its maximum value (rh⁺D) and b_(−j−1) is also maximum (=r−1), a value of q_(−j)=α should guarantee that R_(j)≦h⁺D), thus:

${\left( {{{rh}^{+}D} - {\alpha\; D} + {\frac{r - 1}{r}\omega\; D\mspace{14mu}\ldots}}\mspace{14mu} \right) \leq {h^{+}D}},{then}$ $\alpha \geq {\left( {r - 1} \right){\left( {h^{+} + \frac{\omega}{r}} \right).}}$

Alternatively, we can write

$h^{+} \leq {\frac{\alpha}{r - 1} - {\frac{\omega}{r}.}}$

By replacing ω by ω_(max) in the above equation, we obtain a lower bound expression for h⁺ that guarantees satisfaction of the constraint R_(j)≦h⁺D. Thus, h⁺ is taken as:

$\begin{matrix} {h^{+} = {\frac{\alpha}{r - 1} - {\frac{\omega_{\max}}{r}.}}} & (16) \end{matrix}$

Equation 16 clearly shows that the upper bound for α is (r−1), in which case

$h^{+} = {\left( {1 - \frac{\omega_{\max}}{r}} \right).}$

Likewise, when the shifted partial residue rR_(j-1) equals its minimum value (—rh⁻D) and b_(−j−1) is also minimum (=0), a value of q_(−j)=−α should guarantee that R_(j)≧−h⁻D, thus: (−rh ⁻ D+αD+0)≧−h ⁻ D _(∴)∴α≧(r−1)h ⁻

Alternatively, we can write: h⁻≧α/(r−1)<1.

Thus, to guarantee satisfaction of the constraint R_(j)≧−h⁻D and for maximum overlap regions on the P-D diagram (and accordingly simpler quotient digit selection logic), we use the highest value for h⁻ given by:

$\begin{matrix} {h^{-} = {\frac{\alpha}{\left( {r - 1} \right)} = h}} & (17) \end{matrix}$ where h is the redundancy factor, and we can re-write the equation of h⁺ as follows:

$\begin{matrix} {h^{+} = {{\frac{\alpha}{\left( {r - 1} \right)} - \frac{\omega_{\max}}{r}} = {h - \frac{\omega_{\max}}{r}}}} & (18) \end{matrix}$

In a P-D diagram, we need to determine the selection interval defined by the upper (U_(k)) and lower (L_(k)) bounds of the shifted partial residue (P=rR_(j-1)) for which a given quotient digit value (q_(−j)=k) may be selected such that the next partial residue (R_(j)) satisfies −h⁻D≦R_(j)≦+h⁺D.

From equation (7), we can write P=rR_(j-1)=R_(j)+q_(−j)D−b_(−j−1)Ar⁻¹. Accordingly, we define U_(k) as the upper bound of P (=rR_(j-1)) for which q_(−j)=k yields a valid R_(j) value −h⁻D≦R_(j)≦+h⁺D. Thus:

$\begin{matrix} {{{U_{k} = {{{rR}_{j - 1}\left( {q_{- j} = k} \right)} = {R_{jmax} + {kD} - {b_{{({{- j} - 1})}\max}{Ar}^{- 1}}}}},{or}}{U_{k} = {{h^{+}D} + {kD} - {\left( \frac{r - 1}{r} \right)\omega_{\max}D}}}{U_{k} = {\left( {k + h - \omega_{\max}} \right){D.}}}} & (19) \end{matrix}$

Likewise, we define L_(k) as the lower bound of P (=rR_(j-1)) for which q_(−j)=k yields a valid R_(j) value −h⁻D≦R_(j)≦+h⁺D. Thus: L _(k) =R _(jmin) +kD−b _((−j-1)min) Ar ⁻¹, or L _(k)=(k−h ⁻)D=(k−h)D.  (20)

Using all bits of P and D(2 w+2m+Z+1) bits as input to the quotient digit selection function {SEL(P, D)} requires huge ROM or PLA sizes. For example if w=24 bits, r=8 (i.e., m=3), and Z=6, the quotient digit selection function will have a minimum of 61 input bits, assuming non-redundant representation of both P and D. With such large number of input bits, the hardware complexity of the SEL function is bound to be enormous. For example, if a ROM is used to store this function, the required ROM size (2⁶¹×4-bits) is prohibitively large.

Accordingly, it is advantageous to minimize the number of input bits to the quotient digit selection function. Effectively, we need to use truncated values of P and D, with the smallest possible number of bits as input to the quotient digit selection function. Let these truncated values be P_(t) and D_(t) and let the number of fractional bits of these parameters be n_(p) and n_(D), respectively. Thus, the maximum truncation error values for P and D are 2^(−np) and 2^(−nD) respectively. Using a 2's complement representation, the introduced truncation errors are always positive, i.e., P≧P_(t) and D≧D_(t). We now derive expressions for the optimal values of n_(p), n_(D) and z in terms of the radix r, the redundancy factor h, and the digit set α.

Thus, the selection function defines for each interval of the divisor D [d_(i), d_(i+1)), where d_(i+1)=d_(i)+2^(−nD), comparison constants m_(k)(i) within the overlap region for all values of kε{−(α−1), −(α−2), . . . , −1, 0, +1, . . . , +α}. Since P is represented in the 2's complement system, then P≧P_(t). Accordingly, any given value of P_(t) represents a range of P that is defined by: P_(t)≦P<P_(t)+2^(−nP). Likewise, D≧D_(t). Accordingly, any given value of D_(t) represents a range of D that is defined by D_(t)≦D<D_(t)+2^(−nd).

The set of comparison constants for each range of D is determined such that a given value of P is compared to these constants, based upon which a proper value of q_(−j) is chosen. Thus, for the i^(th) range of D, define the comparison constants m_(k)(i)[k=−(α−1), (α−2), . . . , −1,0,1, . . . , +α] such that: IF m _(k)(i)≦P≦m _(k+1)(i) then q _(−j) =k.  (21)

The P-D diagram is used to help determine these comparison constants. The comparison constants m_(k)(i) are chosen within the overlap regions where a choice of a q_(−j) value of either k or k−1 satisfies the constraint −h⁻D≦R_(j)≦+h⁺D. Since any value within the overlap region may be used as a comparison constant for this region, the choice should be made such that (n_(p)+n_(D)) is minimized.

FIG. 3 shows on a P-D diagram 300 the comparison constants m_(k(i)) for a given truncated divisor value D_(t)=d_(i). As shown in FIG. 3, proper values of the quotient digit for various ranges of P values for the case of r=4 and α=2 can be determined from the P-D diagram 300.

Two conditions must be satisfied when determining the comparison constants: (1) containment, where L_(k)≦m_(k)≦U_(k); and (2) continuity, so that if P=(m_(k)−2^(−nP)), then q_(−j) must equal k−1, which implies that m_(k)−2^(−nP)≦U_(k−1). Written differently, we must have m_(k)≦U_(k−1)+2^(−nP) as well as satisfy the containment constraint. Accordingly, m_(k) should satisfy L_(k)≦m_(k)≦U_(k−1)+2^(−nP).

For a given value of P_(t) the uncertainty in the value of P has an upper bound of ΔP=2^(−nP), i.e., P_(t)≦P<P_(t)+ΔP; accordingly, the upper bound for m_(k) must be reduced by ΔP and accordingly m_(k) should satisfy: L _(k) ≦m _(k) ≦U _(k−1).  (22)

For a feasible m_(k) value, the height of the overlap region (Δy) at a given divisor value (D) must be greater than or equal to the minimum grid 2^(−np); thus Δy=U_(k−1)−L_(k)=(2h−1−ω_(max))D. At D=D_(min), the height of the overlap region Δy is minimum (Δy_(min)), defining the upper bound for 2^(−nP), i.e.,

$\begin{matrix} {{{\Delta\; y_{\min}} = {{U_{k - 1} - L_{k}} = {{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}} > 2^{- {nP}}}}}{or}{2^{nP} \geq {\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}.}}} & (23) \end{matrix}$ For a feasible solution, we must have ω_(max)<(2h−1), i.e., 2^(z)>1/(h−0.5). Thus, the minimum value of n_(p) is given by:

$\begin{matrix} {{{n_{P}\left( \min \right)} = \left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil}{{n_{P}\left( \min \right)} = \left\lceil {{- {{Log}_{2}\left( {{2h} - 1 - \omega_{\max}} \right)}} - {{Log}_{2}\left( D_{\min} \right)}} \right\rceil}{where}{h = {{\alpha/r} - 1.}}} & (24) \end{matrix}$

The lower bound of n_(p)(min) is reached at very high values of Z (Z→∞), in which case (ω_(max)→0) and is given by:

$\begin{matrix} {{n_{P}({Low\_ Bound})} = {\left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1} \right)D_{\min}}} \right\rceil.}} & (25) \end{matrix}$

We define Z₁ as the value of Z at which n_(p)(min) is equal to its lower bound value as follows:

$\begin{matrix} {Z_{1} = {1 + \left\lceil {{Log}_{2}\left\{ \frac{1}{{2h} - 1 - \left( \frac{2^{- {n_{P}{({{Low}\_{Bound}})}}}}{D_{\min}} \right)} \right\}} \right\rceil}} & (26) \end{matrix}$ to determine the overlap region between U_(k−1) and L_(k), based upon which we define the comparison constants that determine the value of the next quotient digit q_(−j) (it should be noted that the negative (P<0) and positive (P>0) overlap regions are not symmetric, and should be considered independently).

The overlap region for a given divisor value, D, is the range of P values where the next quotient digit q_(−j) may be assigned either a value of k−1 or k yielding a value of the next partial residue (R_(j)) which satisfies the range constraint −h⁻D≦R_(j)≦+h⁺D in either case. As defined by equation (22), this range is bounded between U_(k−1) and L_(k). FIG. 4 shows part of the P-D diagram 400 for P>0 and the overlap region where we can select q_(−j)ε{k−1, k}. The comparison constant m_(k)(i) falls in the overlap region 405 between U_(k−1) and L_(k) for the i^(th) divisor range [D:D+2^(−n) ^(D) ). The lower and upper bound values (P_(lower) & P_(upper)) for this comparison constant are given by: P _(lower) =L _(k)(D+2^(−n) ^(D) )=(k−h)(D+2^(−n) ^(D) ) 1≦k≦α  (27) P _(upper) =U _(k−1)(D)=(k+h−1−ω_(max))D. 1≦k≦α  (28) Thus, the selection constants m_(k)(i) are determined for the P>0 range such that m_(k)(i) is an integer multiple of 2^(−nP) and: (k−h)(D+2^(−n) ^(D) )≦m _(k)≦(k+h−1−ω_(max))D 1≦k≦α PositiveOverlap_(k) =Δy ⁺ =P _(upper) −P _(lower)≧0  (29) Δy ⁺=(2h−1−ω_(max))D−(k−h)2^(−n) ^(D) ≧0 1≦k≦α  (30)

For negative overlap regions where P<0, shown in the diagram 500 of FIG. 5, the following equations hold: P _(lower) =L _(k)(D)=(k−h)D −(α−1)≦k≦0  (31) P _(upper) =U _(k−1)(D+2^(−n) ^(D) )=(k+h−1−ω_(max))(D+2^(−n) ^(D) ) −(α−1)≦k≦0  (32) Thus, the selection constants m_(k)(i) are determined for the P<0 range such that m_(k)(i) is an integer multiple of 2^(−nP) and: (k−h ⁻)D≦m _(k)≦(k+h−1−ω_(max))(D+2^(−n) ^(D) ) NegativeOverlap_(k) =Δy ⁻ =P _(upper) −P _(lower)≧0 Δy ⁻=(2h−1−ω_(max))D+(k+h−1−ω_(max))2^(−n) ^(D) ≧0  (33) −(α−1)≦k≦0.  (34)

Valid m_(k(i)) are shown in the shaded regions 505. It should be understood that the present invention contemplates an overlap range (Δy⁺ or Δy⁻) that is smaller for smaller values of D. Higher values of |k| yield smaller overlap regions for both Δy⁺ and Δy⁻. For worst case analysis, the smallest value of D (i.e., D_(min)) and the highest value of |k| (k=α for Δy⁺ or k=−(α−1) for Δy⁻) should be used to yield the minimum overlap. Additionally, worst case analysis is performed on Δy⁻(k=−(α−1), D=D_(min)), since this is where the overlap region is smallest. This yields the following relations:

$\begin{matrix} {2^{n_{D}} \geq \frac{\left( {\alpha - h + \omega_{\max}} \right)}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} & (35) \\ {{n_{D}\left( \min \right)} = {\left\lceil {{Log}_{2}\left\{ \frac{\left( {\alpha - h - \omega_{\max}} \right)}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}} \right\}} \right\rceil.}} & (36) \end{matrix}$

The lower bound of n_(D) is reached at very high values of Z (Z→∞) in which case (ω_(max)→0) and is given by:

$\begin{matrix} {{n_{D}({Low\_ Bound})} = {\left\lceil {{Log}_{2}\left\{ \frac{\left( {\alpha - h} \right)}{\left( {{2h} - 1} \right)D_{\min}} \right\}} \right\rceil.}} & (37) \end{matrix}$ We define Z₂ as the value of Z at which n_(D)(min) is equal to its lower bound value as follows:

$\begin{matrix} {Z_{2} = {1 + {\left\lceil {{Log}_{2}\left\{ \frac{\left( {1 + {D_{\min}2^{n_{d}{({{Low}\_{Bound}})}}}} \right)}{{D_{\min}2^{n_{D}{({{Low}\_{Bound}})}}\left( {{2h} - 1} \right)} - \alpha + h} \right\}} \right\rceil.}}} & (38) \end{matrix}$ The value of Z to be used is the maximum of either Z₁ or Z₂, i.e., Z=MAX(Z ₁ ,Z ₂).  (39)

Carry-Save Adders (CSAs) may be used to evaluate the partial residue (R_(j)) of the recurrence relation, in which case R_(j) is represented in a redundant form as two quantities; a sum and a carry. Assuming 2's complement number representation, and using only n_(p) fractional bits, the truncation error is always positive. A fast carry-propagate adder (CPA) may be used to add the most significant bits of the CSA, which may be used as input to the quotient digit selection function.

Assuming n_(p) to be the number of fractional bits of P used as input to the CPA, the error introduced due to the use of CSA's is less than 2^(−np). In this case, the upper bound for the comparison constants should be reduced by the same amount, and accordingly equation (22) is modified for CSA's to become: L _(k) ≦m _(k) ≦U _(k−1)−2^(n) ^(p) .

The overlap region for a given divisor value, D, is the range of P values where the next quotient digit q_(−j) may be assigned either a value of k−1 or k, yielding a value of the next partial residue (R_(j)) which satisfies the range constraint −h⁻D≦R_(j)≦+h⁺D. FIG. 4 shows part of the P-D diagram for P>0 and the overlap region where we can select q_(−j)ε{k−1, k}. The comparison constant m_(k)(i) falls in the overlap region between (U_(k−1)−2^(−nP)) and L_(k) for a given divisor range [D: D+2^(−n) ^(D) ). The lower and upper bound values (P_(lower) & P_(upper)) for this comparison constant are given by: P _(lower) =L _(k)(D+2^(−n) ^(D) )=(k−h)(D+2^(−n) ^(D) ) 1≦k≦α  (40) P _(upper) =U _(k→1)(D)=(k+h−1−ω_(max))D−2^(−n) ^(P) 1≦k≦α  (41)

On the other hand, when P<0, the bounds are: P _(lower) =L _(k)(D)=(k−h)D −(α−1)≦k≦0  (42) P _(upper) =U _(k−1)(D+2^(−n) ^(D) )=(k+h−1−ω_(max))(D+2^(−n) ^(D) )−2^(−n) ^(P) −(α−1)≦k≦0  (43) Comparing heights of the overlap regions for CSA's (equations 40-43) and CPA's (equations 31-34), it can be seen that the overlap region height for CSA's is lower by 2^(−n) ^(p) due to the lowered upper boundary of this region (P_(upper)) by the same amount.

With CSA's, the selection constants m_(k)(i) are determined such that m_(k)(i) is an integer multiple of 2^(−n) ^(p) and: (k−h)(D+2^(−n) ^(D) )≦m _(k)≦(k+h−1−ω_(max))D−2^(−n) ^(p) 1≦k≦α (P>0) (k−h)D≦m _(k)≦(k+h−1−ω_(max))(D+2^(−n) ^(D) )−2^(n) ^(P) −(α−1)≦k≦0 (P<0)

To derive mathematical expressions for optimal parameter values, we consider the worst case overlap region, i.e., Δy_(min) ⁻, at the smallest value of D (i.e., D_(min)) and the most negative value of k (i.e., k=−(α−1)).

$\begin{matrix} {{\Delta\; y_{\min}^{-}} = {{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}} > {2^{- {nP}} - {\left( {\alpha - h + \omega_{\max}} \right)2^{- n_{D}}}} \geq 0}} & (44) \\ {\mspace{79mu}{2^{n_{D}} \geq \frac{a - h + \omega_{\max}}{{\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}} - 2^{- n_{P}}}}} & (45) \end{matrix}$ To have a feasible solution, the denominator of equation (45) must be greater than zero. Thus:

$\begin{matrix} {{{2^{- n_{P}} < {\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}}},{and}}{{n_{P}\left( \min \right)} = \left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil}} & (46) \end{matrix}$ Thus, the minimum feasible value of n_(P) is given by:

$\begin{matrix} {{n_{P}\left( \min \right)} = \left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - 2^{1 - z}} \right)D_{\min}}} \right\rceil} & (47) \end{matrix}$

The lower bound of the minimum n_(P) is reached at very high values of z (z→∞) and is given by:

$\begin{matrix} {{n_{P}({Low\_ Bound})} = {\left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1} \right)D_{\min}}} \right\rceil.}} & (48) \end{matrix}$ Let z₁ be the minimum value of z at which n_(p(min))=n_(P)(Low_Bound), thus:

$\begin{matrix} {z_{1} = {\left\lceil {1 - {{Log}_{2}\left\{ {{2h} - 1 - \frac{2^{- {n_{P}{({{Low}\_{Bound}})}}}}{D_{\min}}} \right\}}} \right\rceil.}} & (49) \end{matrix}$ Multiplying both sides of equation (45) by 2^(n) ^(P)

$\begin{matrix} {2^{n_{P} + n_{D}} \geq {\frac{\alpha - h + \omega_{\max}}{{{D_{\min}\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack}2^{- n_{P}}} - 2^{{- 2}n_{P}}}.}} & (50) \end{matrix}$

To reduce the complexity of the quotient digit selection hardware, (n_(P)+n_(D)) must be minimized or, equivalently, the value of 2^(n) _(p) ^(+n) _(D) must be minimized. Defining Y=2^(n) ^(p) ^(+n) ^(D) and X=2^(−n) ^(p) , equation (50) is rewritten as

$\begin{matrix} {Y \geq {\frac{\alpha - h - \omega_{\max}}{{{D_{\min}\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack}X} - X^{2}}.}} & (51) \end{matrix}$ Differentiating equation (51) with respect to X and equating

${\frac{\partial Y}{\partial X} = 0},$ we get:

$\begin{matrix} {{X = {2^{- n_{P}} = {{\frac{1}{2}\left\lbrack {\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}} \right\rbrack} = {\left( {h - 0.5 - 2^{- t}} \right){D_{\min}.{Thus}}}}}},} & (52) \\ {{n_{p}({opt})} = {1 + {\left\lceil {{Log}_{2}\frac{1}{\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}}} \right\rceil.}}} & (53) \end{matrix}$

In general,

${Log}_{2}\frac{1}{\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}}$ yields a non-integer value. Accordingly, the actual optimal n_(p) value may be either the rounded up or rounded down to the integer value nearest to the value computed by equation (53). Since equation (53) yields a value for n_(p) that is higher only by 1 than the minimum n_(p) value defined by equation (46), it is clear that the optimal n_(p) value may either equal the minimum value specified by equation (46), or it may be larger by just one bit. The optimum no value may be computed from equation (45):

$\begin{matrix} {\mspace{79mu}{{2^{n_{D}} \geq \frac{2\left( {\alpha - h + \omega_{\max}} \right)}{\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}}}{{n_{D}({opt})} = {1 + \left\lceil {{{Log}_{2}\left\{ {\alpha - h + \omega_{\max}} \right\}} - {{Log}_{2}\left\{ {\left\lbrack {{2h} - 1 - \omega_{\max}} \right\rbrack D_{\min}} \right\}}} \right\rceil}}}} & (54) \end{matrix}$

The lower bound of the optimal n_(D) value is reached at very high values of z (z→∞) and is given by: n _(D)(Low_Bound)=1+|Log₂ {α−h}−Log₂{(2h−1)D _(min)}|.  (55)

Let z₂ be the minimum value of z at which n_(D)(opt)=n_(D)(Low_Bound), Using equation (54), z₂ can be derived as follows:

$\begin{matrix} {z_{2} = {1 + {\left\lceil {{Log}_{2}\left\{ \frac{{2^{{n_{D}{({Low\_ Bound})}} - 1}D_{\min}} + 1}{{2^{{n_{D}{({Low\_ Bound})}} - 1}{D_{\min}\left\lbrack {{2h} - 1} \right\rbrack}} - \alpha + h} \right\}} \right\rceil.}}} & (56) \end{matrix}$ The optimal value of z(z_(opt)) is the larger of z₁ or z₂; thus: Z _(opt)=MAX(z ₁ ,z ₂).  (57)

The value of Z is chosen as the higher of two values, Z₁ and Z₂ that are derived from the lower bound values of n_(p) and n_(D). Expressions for Z₁ and Z₂ have been derived for the case where carry-propagate adders are used, as well as the case where carry-save adders are used.

Whether carry-propagate or carry-save adders are used, the value of Z₁ is the same, since the expressions for n_(P)(Low_Bound) and Z₁ are identical in both cases, as can be readily seen by comparing equations (25) and (26) on the one hand with equations (48) and (49) on the other.

For the case of carry-save adders, the low bound value of n_(D) as given by equation (54) is higher by 1 than its value for the carry-propagate adder case as given by equation (37). Accordingly, the value of Z₂ for both the CPA and CSA cases as computed by equations (38) and (55), respectively, are the same.

Accordingly, the value of Z is independent of the type of adder that is used for implementation. Thus, the equations to derive Z are summarized below:

${n_{P}({Low\_ Bound})} = \left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1} \right)D_{\min}}} \right\rceil$ $z_{1} = {1 + \left\lceil {{Log}_{2}\left\{ \frac{1}{\left\{ {{2h} - 1 - \left( \frac{2^{- {n_{P}{({low\_ Bound})}}}}{D_{\min}} \right)} \right\}} \right\}} \right\rceil}$ ${n_{D}({Low\_ Bound})} = \left\lceil {{Log}_{2}\left\{ \frac{\left( {\alpha - h} \right)}{\left( {{2h} - 1} \right)D_{\min}} \right\}} \right\rceil$ $Z_{2} = {1 + \left\lceil {{Log}_{2}\left\{ \frac{\left( {1 + {D_{\min}2^{n_{D}{({Low\_ Bound})}}}} \right)}{{D_{\min}2^{n_{D}{({Low\_ Bound})}}\left( {{2h} - 1} \right)} - \alpha + h} \right\}} \right\rceil}$ Z = MAX(Z₁, Z₂).

It should be noted that at α=r−1, the redundancy factor h=1, and the equation of Z₁ yields an infinite value. Likewise, the expression of Z₂ may yield infeasible values for certain cases, e.g., r=4 and α=2. Table II, Table III, and Table IV show the values of Z₁ and Z₂ for several radixes r=2^(m) at α=r/2 (minimal redundancy), α=r−2, and α=r−1 (maximal redundancy).

TABLE II Values of Z₁ and Z₂ for minimal redundancy (α = r/2) m 2 3 4 5 6 r = 2^(m) 4 8 16 32 64 α = r/2 2 4 8 16 32 Z₁ 5 7 9 11 13 Z₂ X 6 8 10 12

TABLE III Values of Z₁ and Z₂ for α = r − 2 m 2 3 4 5 6 r = 2^(m) 4 8 16 32 64 α = r − 2 2 6 14 30 62 Z₁ 5 4 3 3 3 Z₂ X 5 6 7 8

TABLE IV Values of Z₁ and Z₂ for maximal redundancy (α = r − 1) m 2 3 4 5 6 r = 2^(m) 4 8 16 32 64 α = r − I 3 7 15 31 63 Z₁ X X X X X Z₂ X 4  5  6  7

The tables show that the maximum Z value occurs under minimal redundancy conditions (α=r/2), and is equal to (2 m+1). For cases where the expression of either Z₁ or Z₂ yields an infeasible value, some choice criterion may be used to define the value of Z. For example, Z may be chosen equal to the feasible value of either Z₁ or Z₂ totally neglecting the one with infeasible value. An alternative strategy is: if Z₁ has an infeasible value, increment the value of n_(P)(Low_Bound) by 1 and recompute Z₁; or if Z₂ has an infeasible value, increment the value of n_(D)(Low_Bound) by 1 and recompute Z₂.

Alternative approaches are also possible, e.g., adopting values of the closest higher system, e.g., for the case of Z₂ with r=4 and α=r−2=2, we set Z₂=5 corresponding to the system with r=8 and α=r−2=6.

Based on the above developed theory, given the system radix r, and the quotient digit set parameter α, the optimal parameters for the high-radix multiplier-divider may be determined as follows:

${1\text{-}{Compute}\mspace{14mu} h} = \frac{\alpha}{r - 1}$

-   -   2—Determine Z as detailed supra     -   3—Compute ω_(max)=2^(1-z)

${4\text{-}{Compute}\mspace{14mu}{n_{P}\left( \min \right)}} = \left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil$

-   -   5—For n_(p)=n_(p)(min) To n_(p)(min)+2, repeat:         -   a. Compute n_(D)(min) corresponding to the current n_(p)             value (equation 45 for CSA, or equation 35 for CPA)         -   b. For the current n_(p) and n_(D) values, repeat the             following until either a solution is found or for a maximum             of 3 times;             -   i. Determine the comparison constants M_(k)(i)*             -   ii. if a feasible solution exists (i.e., a complete set                 of comparison constants is obtained) then store                 n_(p)+n_(D), else increment n_(D)     -   6—Select the solution which yields the smallest n_(p)+n_(D). In         case of more than one solution having the same smallest         n_(p)+n_(D), choose the one with the smallest n_(p).

The comparison constants m_(k)(i) are determined to satisfy the following:

-   -   1. m_(k)(i) is an integer multiple of 2^(−nP);     -   2. L_(k) (D+2^(−n) ^(D) )≦m_(k)≦U_(k−1)(D) 1≦k≦α (P>0)     -   3. L_(k) (D)≦m_(k)≦U_(k−1)(D+2^(−n) ^(D) ) −(α−1)≦k≦0 (P<0)     -   where L_(k)=(k−h)D, and         -   U_(k−1)(k+h−1−ω_(max))D−2^(−n) ^(p) if Carry-Save Adders are             used; or         -   =(k+h−1−ω_(max)) D if Carry-Propagate Adders are used.

FIG. 1 shows an embodiment of the high-radix multiplier-divider 100 that utilizes carry-save adders. High-radix multiplier-divider 100 has the following features. A counter 101 is used to hold the number of iterations to be performed

$\left( {\eta_{iterate} = \left\lceil \frac{k + Z}{m} \right\rceil} \right).$ The m most significant bits of the k-bit B-register constitute the current digit (b_(−j−1)) of B. In each iteration register B is shifted left by m-bits. The selection function is implemented either as a ROM 105 or a PLA where the truncated values of P and D (i.e., P_(t) and D_(t)) are the input to this ROM 105 (or PLA) for a total of (n_(D)+n_(P)+m) bits. The output of the ROM/PLA 105 is the (m+1)-bit signed value of q_(−j).

The value of P(=rR_(j)) uses a redundant representation in the form of a SUM component (PS), and a CARRY component (PC), which are held in the registers PSR 120 and PCR 118, respectively. Accordingly there are four quantities that need to be added in each iteration (i.e., each execution of the recurrence relation of equation (7) R_(j)=rR_(j-1)−q_(−j+1)D+b_(−j)Ar⁻¹; namely PS, PC, (−q_(−j)*D), and (b_(−j−1)*A/r).

The multiplexer MUXa 110 generates the k+m bits (b_(−j−1)*A′), which is left appended by 1+Z+m bits of 0 value to generate the signed quantity (b_(−j−1)*A/r), where A=A′/2^(z). The multiplexer MUXd 107 generates k+m+1 bits of the signed 1's complement of the quantity (−q_(−j)*D), i.e., if q_(−j) is positive, ( q_(−j)D) is generated, while if it is negative, (|q_(−j)|*D) is generated. The output of MUXd 107 is left appended by Z+m bits, each having a value that equals the sign bit of (q_(−j)D) (sign-bit extension).

A Carry-Lookahead Adder 135 (CLA) is used to add the (1+m+n_(P)) most significant bits of the sum and carry components of the shifted partial residue (PS & PC). The resulting summation is the truncated P_(t) value used as input to the ROM/PLA 105. Adding the 4 quantities PS, PC, (−q_(−j)*D), and (b_(−j−1)*A/r) is done using two Carry-Save adders, CSA1 112 and CSA2 115. CSA1 112 adds PS, PC, and (b_(−j−1)*A/r) yielding two outputs: a partial sum component (Sum 1), and a partial carry component (Cry 1). The second CSA, CSA2 115 adds Sum1, Cry1, and (−q_(−j)*D). For a correct result, the 1's complement representation of (−q_(−j)*D) is turned into a 2's complement by forcing bit 0 of Cry 1 to equal the sign bit value of (−q_(−j)*D). CSA2 115 yields two outputs; a partial sum component (Sum2), and a partial carry component (Cry2).

An m-bit left-shifted version of Sum2 and an m-bit left-shifted version of Cry2 are stored in two registers, PSR 120 and PCR 118 to represent (rR_(j)). The outputs of PSR 120 and PCR 118 are fed-back as input to CSA1 112, representing the shifted partial residue (rR_(j)), while the (1+m+n_(P)) most significant bits of PSR 120 and PCR 118 are added using the CLA 135 to yield the value of P_(t).

At the last iteration, a second CLA 125 is used to assimilate the sum and carry components of the shifted partial residue (PS & PC) to yield the value of P. This CLA 125 may via last cycle AND gate 130, or, alternatively, may not utilize the (1+m+n_(P))-bit first CLA 135 as part of it to yield the (1+m+n_(P)) most-significant bits of the result, as shown in FIG. 1.

Example 1

If Carry-Propagate Adders (CPA) are used, compute AB/D=(0.1110_(—)1001)*(0.10011100)/(0.1101_(—)1110) using radix r=4, and D_(min)=0.5. Let A′=0.1110_(—)1001, B=0.10011100, and D=0.1101_(—)1110.

h=⅔=0.667, n_(P)(Low_Bound)=3, n_(D)(Low_Bound)=3, Z₁=5, Z₂ computed at [nD(Low_Bound)+1] is 4. Thus, Z=Max(Z₁, Z₂)=5, ω_(max)=2⁻⁴=0.0625,

${n_{P}\left( \min \right)} = {\left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil = 3}$ and ${n_{D}\left( \min \right)} = {\left\lceil {{Log}_{2}\left\{ \frac{\left( {\alpha - h + \omega_{\max}} \right)}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}} \right\}} \right\rceil = 4.}$

Considering the worst case of D=D_(min)=0.5, and computing L_(k) and U_(k−1) at various values of k, it can be shown that no solution is possible for n_(P)=3 and n_(D)=4. However, a solution exists for the case of n_(P)=3 and n_(D)=5. The table below lists the values of L_(k) and U_(k−1) for various values of k at D=D_(min)=0.5, in addition to possible values of the comparison constants m_(k) for this case.

TABLE V L_(k), U_(k−1), and m_(k) for D = D_(min) = 0.5 k L_(k) U_(k−1) m_(k) −1 −0.83333 −0.74154 −0.75 0 −0.33333 −0.21029 −0.25 1 0.177083 0.302083 0.25 2 0.708333 0.802083 0.75

For the example at hand since n_(D)=5, the truncated value of D is given by D_(t)=0.11011=27/32. Table VI gives the computed values of L_(k) and U_(k−1) for various values of k, together with the selected values of comparison constants for the range D=27/32:28/32.

TABLE VI L_(k), U_(k−1) and m_(k) for D = 27/32 k L_(k) U_(k−1) m_(k) −1 −1.40625 −1.22135 −1.25 0 −0.5625 −0.34635 −0.375 1 0.291667 0.509766 0.375 2 1.166667 1.353516 1.25

The pre-processing steps are given below:

TABLE VII Preliminary values A = A′ * 2^(−z) 0.0_0000_1110_1001 A/r 000.0000_0001_1101 001 A * B = 2⁻⁵ * (0.1000 1101 1111 11)

In Tables VIII and IX, below, values of b_(−j) are noted as b(−1), b(−2) etc., values of q_(−j) are notes as q(−1), q(−2), etc., and values of R_(j) are noted as R0, R1, etc.

TABLE VIII Values of b_(−j) b(−1) = 2 b(−2) 1 b(−3) 3 b(−4) 0 b(−5:−8) 0

TABLE IX Calculation of R₁-R₇ and q⁻¹-q⁻⁷ Processor Size 8 + (2*2) + 5 + 1 = 18 bits # iterations = 4 + [5/2] = 7 R0 = b(−1)A/r 000.0000_0011_1010_010 rR0 000.0000_1110_1001_000 → q(−1) = 0 +b(−2)*A/r 000.0000 0001 1101 001 R1 000.0001 _0000_01 1 0_001 rR1 000.0100_0001_1000_100 → q(−2) = 0 +b(−3)*A/r 000.0000 0101 0111 011 R2 000.0100_0110_1111_111 rR2 001.0001_1011_1111_100 → q(−3) = 1 −q(−3)D 111.0010_0010_0000_000 +b(−4)*A/r 000.0000 0000 0000 000 R3 000.0011_1101_1111_100 rR3 000.1111_0111_1110_000 → q(−4) = 1 −q(−4)D 111.0010_0010_0000_000 +b(−5)*A/r 000.0000 0000 0000 000 R4 000.0001_1001_1110_000 rR4 000.0110_0111_1000_000 → q(−5) = 1 ·q(−5)D 111.0010_0010_0000_000 +b(−6)*A/r 000.0000 0000 0000 000 R5 111.1000_1001_1000_000 rR5 110.0010_0110_0000_000 → q(−6) = 2 −q(−6)D 001.1011_1100_0000_000 +b(−7)*A/r 000.0000 0000 0000 000 R6 111.1110_0010_0000_000 rR6 111.1000_1000_0000_000 → q(−7) = −1 −q(−7)D 000.1101_1110_0000_000 +b(−8)*A/r 000.0000 0000 0000 000 R7 000.0110_0110_0000_000

From the above, the resulting quotient may be expressed as: Q=[0.00111(−2)(−1)]₄=0.0000_(—)0_(—)1010_(—)0011_(—)1=2⁻⁵*(0.1010_(—)0011_(—)1)₂ In the foregoing, it is noted that Q has nine significant bits, whereas only eight bits are required. The remainder may be expressed as: R ₇=000.0110_(—)0110_(—)0000_(—)000 R=r ⁻⁷ R ₇=2⁻¹⁵*(0.1100_(—)1100)

The validity of the result may be verified from the following considerations: A*B=2⁻⁵*(0.1000_(—)1101_(—)1111_(—)11) Q*D=2⁻⁵*(0.1000_(—)1101_(—)1100_(—)1001) R=r ⁻⁷ R ₇=2⁻⁵*2⁻¹⁰*(0.1100_(—)1100) Q*D+R=A*B=2⁻⁵*(0.1000_(—)1101_(—)1111_(—)11) Since the required accuracy is only 8 bits, a correction step is required. In the correction step, Q=Q−ulp=2⁻⁵*(0.1010_(—)0011)₂. But, AB=DQ+R=D*(Q−ulp)+(R+D*ulp). Therefore, the corrected quotient is Q=Q−ulp=2⁻⁵*(0.1010_(—)0011) and the corrected remainder is: R=R+D*ulp=2⁻¹⁵*(0.1100_(—)1100)+2⁻¹⁴*(0.1101_(—)1110)=2⁻¹³*(0.1010_(—)0010)

In a required correction step, since the original operand is A′=0.1110_(—)1001=A*2^(z)=A*2⁵, the actual quotient Q′ and remainder R′ are given by: Q′=Q*2⁵=0.1010_(—)0011 R′=R*2⁵=2⁻⁸*(0.1010_(—)0010)

Example 2

Using Carry-Save Adders (CSA), compute AB/D=(0.1110_(—)1001)*(0.10011100)/(0.1101_(—)1110) using radix r=4, D_(min)=0.5.

Let A′=0.1110_(—)1001, B=0.10011100, and D=0.1101_(—)1110. h=⅔=0.667, n_(P)(Low_Bound)=3, n_(D)(Low_Bound)=4, Z₁=5, Z₂=4. Thus,

${Z = {{{Max}\left( {{Z\; 1},{Z\; 2}} \right)} = 5}},{\omega_{\max} = {2^{- 4} = {.0625}}},{{n_{P}\left( \min \right)} = {\left\lceil {{Log}_{2}\frac{1}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil = 3}},{and}$ ${n_{D}\left( \min \right)} = {\left\lceil {{Log}_{2}\frac{\left( {\alpha - h + \omega_{\max}} \right)}{\left( {{2h} - 1 - \omega_{\max}} \right)D_{\min}}} \right\rceil = 4.}$

Considering the worst case of D=D_(min)=0.5, and computing L_(k) and U_(k−1) at various values of k, it can be shown that no solution is possible for n_(p)=3 and n_(D)=4. However, a solution exists for the case of n_(p)=4 and n_(D)=6. Table X below lists the values of L_(k) and U_(k−1) for various values of k at D=D_(min)=0.5, in addition to possible values of the comparison constants m_(k) for this case.

TABLE X L_(k), U_(k−1), and m_(k) for n_(P) = 4 and n_(D) = 6 k L_(k) U_(k−1) m_(k)(0.5) −1 −0.83333 −0.78223 −0.8125 0 −0.33333 −0.2666 −0.3125 1 0.171875 0.239583 0.1875 2 0.6875 0.739583 0.6875

For the example at hand since n_(D)=6, the truncated value of D is given by D_(t)=0.110111=55/64. Table XI gives the computed values of L_(k) and U_(k−1) for various values of k, together with the selected values of comparison constants for the range D=55/64:56/64.

TABLE XI L_(k), U_(k−1), and m_(k) for D = 55/64 k L_(k) U_(k−1) m_(k)(55/64) −1 −1.43229 −1.28385 −1.3125 0 −0.57292 0.40885 −0.4375 1 0.291667 0.456706 +0.4375 2 1.166667 1.316081 1.3125

Interim calculated values for each iteration are shown in the table 700 in FIG. 7. Using CPA to add PSR and PCR (as shown in FIG. 1), we obtain R₇=000.0110_(—)0110_(—)0000_(—)000. The resulting quotient, Q, may be expressed as: Q=[0.001102(−1)]₄=0.0000_(—)0_(—)1010_(—)0011_(—)1=2⁻⁵*(0.1010_(—)0011_(—)1)₂ In the foregoing, it is noted that Q has nine significant bits, whereas only eight bits are required. The remainder may be expressed as: R ₇=000.0110_(—)0110_(—)0000_(—)000 R=r ⁻⁷ R ₇=2⁻¹⁵*(0.1100_(—)1100)

The validity of the result may be verified from the following considerations: A*B=2⁻⁵*(0.1000_(—)1101_(—)1111_(—)11) Q*D=2⁻⁵*(0.1000_(—)1101_(—)1100_(—)1001) r=r ⁻⁷ R ₇=2⁻⁵*2⁻¹⁰*(0.1100_(—)1100) Q*D+R=A*B=2⁻⁵*(0.1000_(—)1101_(—)1111_(—)11) The required accuracy is 8 bits, but the computed accuracy is 9 bits. Therefore, the following correction step is needed: Q=Q−ulp=2⁻⁵*(0.1010_(—)0011)₂. But, AB=DQ+R=D*(Q−ulp)+(R+D*ulp). Therefore, the corrected quotient is Q=Q−ulp=2⁻⁵*(0.1010_(—)0011) and the corrected remainder is: R=R+D*ulp=2⁻¹⁵*(0.1100_(—)1100)+2⁻¹⁴*(0.1101_(—)1110)=2⁻¹³*(0.1010_(—)0010)

In a required correction step, since the original operand is A′=0.1110_(—)1001=A*2^(z)=A*2⁵, the actual quotient Q′ and remainder R′ are given by: Q′=Q*2⁵=0.1010_(—)0011 R′=R*2⁵=2⁻⁸*(0.1010_(—)0010)

Example 3

Compute AB/D=(1.110_(—)1001)*(1.001_(—)1100)/(1.101_(—)1110) using radix r=4, D_(min)=1.0 and Carry-Save Adders.

Let A′=1.110_(—)1001, B=1.001_(—)1100 and

$D = {{1.101\_ 1110} = {{1\frac{94}{128}} = {1.734375.}}}$ h=⅔=0.667, n_(P)(Low_Bound)=2, n_(D)(Low_Bound)=3, Z₁=5, Z₂=4, Z=Max(Z₁, Z₂)=5, ω_(max)=2⁻⁴=0.0625.

Considering the worst case of D=D_(min)=1.0, and computing L_(k) and U_(k−1) at various values of k, it can be shown that a solution exists for the case of n_(p)=3 and n_(D)=5. Table XII below lists the values of L_(k) and U_(k−1) for various values of k at D=D_(min)=1.0, in addition to possible values of the comparison constants m_(k) for this case.

TABLE XII L_(k), U_(k−1), and m_(k) for D = D_(min) = 1.0 k L_(k) U_(k−1) m_(k) −1 −1.66667 −1.56445 −1.625 0 −0.66667 −0.5332 −0.625 1 0.34375 0.479167 0.375 2 1.375 1.479167 1.375

For the example at hand, since n_(D)=5, the truncated value of D is given by D_(t)=1.1011=1 11/16. Table XIII gives the computed values of L_(k) and U_(k−1) for various values of k, together with the selected values of comparison constants for the range

$D = {{1\frac{11}{16}}:{1{\frac{12}{16}.}}}$

TABLE XIII L_(k), U_(k−1), and m_(k) for D = 1 11/16 k L_(k) U_(k−1) m_(k)(27/32) −1 −2.8125 −2.52409 −2.625 0 −1.125 −0.80534 −1.0 1 0.572917 0.894531 0.875 2 2.291667 2.582031 2.5

Processor size is 8+(2*2)+5+1=18 bits. The number of iterations is 4+┌5/2┐=7. The results of intermediate calculations per iteration are shown in the table 800 of FIG. 8.

Compared to a pure divider, the size requirements of the quotient digit selection logic may be larger in the multiplier-divider 100 due to the reduced overlap regions in its P-D diagrams. Overall, it is expected that the area of multiplier-divider 100 will be slightly larger than that of a divider only.

It should be noted that the circuit of FIG. 1 is a partial circuit and would also include, e.g., a register for storing the quotient digits selected from the lookup table 105. It will further be noted that the adder 125 of FIG. 1 may be a carry-propagate adder.

As noted above, the circuit may be incorporated into the architecture of a computer processor, into a security coprocessor integrated on a motherboard with a main microprocessor, into a digital signal processor, into an application specific integrated circuit (ASIC), or other circuitry associated with a computer, electronic calculator, or the like. It should be understood that the calculations may be performed by any suitable computer system, such as that diagrammatically shown in FIG. 9. Data is entered into system 100 via any suitable type of user interface 116, and may be stored in memory 112, which may be any suitable type of computer readable and programmable memory. Calculations are performed by processor 114, which may be any suitable type of computer processor and may be displayed to the user on display 118, which may be any suitable type of computer display.

Processor 114 may be associated with, or incorporated into, any suitable type of computing device, for example, a personal computer or a programmable logic controller. The display 118, the processor 114, the memory 112 and any associated computer readable recording media are in communication with one another by any suitable type of data bus, as is well known in the art.

Examples of computer-readable recording media include a magnetic recording apparatus, an optical disk, a magneto-optical disk, and/or a semiconductor memory (for example, RAM, ROM, etc.). Examples of magnetic recording apparatus that may be used in addition to memory 112, or in place of memory 112, include a hard disk device (HDD), a flexible disk (ED), and a magnetic tape (MT). Examples of the optical disk include a DVD (Digital Versatile Disc), a DVD-RAM, a CD-ROM (Compact Disc-Read Only Memory), and a CD-R (Recordable)/RW.

It should be noted that, in the above method, the multiplicand is both arbitrary and input to the system.

It is to be understood that the present invention is not limited to the embodiment described above, but encompasses any and all embodiments within the scope of the following claims. 

We claim:
 1. A high-radix multiplier-divider for performing simultaneous multiplication and division in radix r of a multiplicand A=0·a⁻¹a⁻² . . . a_(−n), a multiplier B=0·b⁻¹b⁻² . . . b_(−n), and a divisor D=0·d⁻¹d⁻² . . . d_(−n) so that $\left( \frac{AB}{D} \right)$ where AB<D, comprising: a plurality of registers for storing the multiplicand, the multiplier, the divisor, a quotient, and a remainder, wherein the multiplicand is an input value; a first data switch having a first input receiving the multiplicand right-shifted by one digit (m-bits), a second input for sequentially receiving bits from the multiplier register corresponding to digits of the multiplier from the most significant digit to the least significant digit, and an output for outputting b_(−j−1)A_(r) ⁻¹ where j is a counter of the multiplier digits; a digital lookup table having a first input for receiving a truncated shifted partial remainder, a second input for receiving a truncated divisor, and an output for outputting a quotient digit corresponding to the truncated shifted partial remainder and truncated partial divisor, the output being stored in the quotient register and left-shifted by the radix; a second data switch having a first input connected to the output of the digital lookup table, a second input receiving the divisor, and an output for outputting the product of the quotient digit and the divisor; and an addition module having inputs for receiving the output of the first data switch and the 2's complement output of the second data switch, the addition module being configured for recursively adding the inputs for each digit of the multiplier and outputting the truncated partial remainder at each iteration for input to the digital lookup table.
 2. The high-radix multiplier-divider according to claim 1, wherein said addition module comprises at least two carry-propagate addition circuits.
 3. The high-radix multiplier-divider according to claim 1, wherein said addition module comprises at least two cascaded carry-save adder circuits, at least one carry-lookahead adder circuit for computing the truncated partial remainder, and at least one carry-propagate adder for adding partial sum bits and partial carry bits from the cascaded carry-save adders on the final iteration, the carry-propagate adder having an output connected to the remainder register.
 4. The high-radix multiplier-divider according to claim 1, wherein said digital lookup table constrains the quotient digit so that the absolute value of the partial remainder is less than the absolute value of the divisor.
 5. The high-radix multiplier-divider according to claim 1, wherein said digital lookup table comprises an area of ROM (read-only memory) having the table stored therein.
 6. The high-radix multiplier-divider according to claim 1, wherein said digital lookup table comprises a programmable logic array.
 7. A computer software product that includes a non-transitory storage medium readable by a processor, the non-transitory storage medium having stored thereon a set of instructions for performing high-radix multiplication and division of fractions in radix r including a multiplicand A=0·a⁻¹a⁻² . . . a_(−n), a multiplier B=0·b⁻¹b⁻² . . . b_(n), and a divisor D=0·d⁻¹d⁻² . . . d_(−n) so that $\left( \frac{AB}{D} \right)$ where AB<D, the instructions comprising: (a) a first set of instructions which, when loaded into main memory and executed by the processor, causes the processor to receive an input value of the multiplicand and initialize a partial remainder to equal b⁻¹Ar⁻¹; (b) a second set of instructions which, when loaded into main memory and executed by the processor, causes the processor to initialize an iteration counter j=1; (c) a third set of instructions which, when loaded into main memory and executed by the processor, causes the processor to initialize a quotient register to equal 0; (d) a fourth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to shift the quotient register left by the radix; (e) a fifth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to look up a quotient digit in a digital lookup table using the partial remainder and the divisor D; (f) a sixth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to store the quotient digit in the quotient register's least significant bits; (g) a seventh set of instructions which, when loaded into main memory and executed by the processor, causes the processor to shift the partial remainder left by the radix; (h) an eighth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to set the partial remainder equal to the shifted partial remainder minus the product of the quotient digit and the divisor plus the quantity b_(−j−1)A_(r) ⁻¹; (i) a ninth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to increment the counter j by one; (j) a tenth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to repeat the fourth set of instructions through the ninth set of instructions by a number of iterations equal to the number of bits in the quotient register divided by Log₂ r; (k) an eleventh set of instructions which, when loaded into main memory and executed by the processor, causes the processor to subtract one from the quotient register and adding the divisor to the partial remainder when the partial remainder of the eighth set of instructions is negative; and (l) a twelfth set of instructions which, when loaded into main memory and executed by the processor, causes the processor to set the quotient equal to the quotient register and the remainder equal to the partial remainder. 